Vortex-lattice melting in two-dimensional superconductors in intermediate fields 
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To examine the field dependence of the vortex lattice melting transition in two-dimensional (2D) 
superconductors, Monte Carlo simulations of the 2D Ginzburg-Landau (GL) model are performed 
by extending the conventional lowest Landau level (LL) approximation to include several higher LL 
modes of the superconducting order parameter with LL indices up to six. It is found that a nearly 
vertical melting line in lower fields, which is familiar within the elastic theory, is reached just by 
including higher LL modes with LL indices less than five, and that the first order character of the 
melting transition in higher fields is significantly weakened with decreasing the field. Nevertheless, 
a genuine crossover to the consecutive continuous melting picture intervened by a hexatic liquid is 
not found within the use of the GL model. 

PACS numbers: 



I. INTRODUCTION 



The vortex phase diagram in type II superconductors 
has been extensively studied in relation to the high T c 
cuprates in magnetic fields which is a typical three di- 
mensional (3D) system with strong fluctuation. The vor- 
tex lattice melting transition in 3D systems in clean limit 
has been examined as a first step for understanding phe- 
nomena in real superconductors with quenched disorder 
and is now believed to be of first order in any magnetic 
field. In contrast, understanding of the phase diagram in 
2D case in nonzero magnetic fields have not progressed 
sufficiently This is partly because the resulting vortex 
lattice itself is nonsuperconducting. No state with zero 
resistance is reached [l[ in clean limit. Further, a weak 
but nonvanishing disorder destroys the quasi long range 
order of vortex positions, and a vortex glass, i.e, a su- 
perconducting vortex phase, is never realized at finite 
temperatures, implying that there will be no phase tran- 
sition in real 2D superconductors with disorder in finite 
fields and finite temperatures 0, H| . 

However, the field-temperature vortex phase diagram 
in a 2D superconductor remains unresolved even theoret- 
ically. Based on the elastic theory, there are at least two 
possibilities, a direct first order melting and the defect- 
unbinding melting composed of two consecutive continu- 
ous transitions and a hexatic liquid crystal phase inter- 
vening between them In high fields where the pair- 
field -0 is limited to the modes in the lowest Landau level 
(LL), however, the 2D melting transition is of first or- 
der according to the direct Monte Carlo simulation of 
the Ginzburg-Landau (GL) model 0, Q . Then, it will be 
valuable to clarify whether this first order transition is 
changed to the consecutive continuous melting scenario 
in lower fields or not. In fact, the character of the melt- 
ing transition has been assumed in previous studies to 
be independent of the strength H of the magnetic field. 
In addition, when considering the real disordered case, 
clarifying this issue on the field dependence of the melt- 



ing mechanism would improve understanding of the 2D 
vortex states in the following sense: In real superconduct- 
ing thin films with weak quenched disorder, a first order 
melting in clean limit is not suggested in thermodynamic 
and resistive data. Note that such experiments are usu- 
ally performed in much lower fields than H C 2(0) where 
roles of quenched disorder are believed to be weaker 0]. 
The fact that no first order melting has been suggested 
so far in physical quantities in real 2D thin films might 
be understood if the melting transition in clean limit is 
continuous or a highly weak first order one in such lower 
fields. 

Based on such a background on the 2D vortex states, 
in the present work, the previous numerical analysis lim- 
ited to the LLL modes of ip is extended to the case with 
several higher LL modes to address the nature of the 
vortex lattice melting in the GL model in lower fields. 
This paper is organized as follows. In sec. II. the model 
and the procedures used for the analysis and simulations 
are explained. In sec. Ill, the obtained numerical results 
are explained and discussed. Section IV includes a sum- 
mary of the present work and a comparison with other 
theoretical works. 



II. MODEL AND PROCEDURES 

Our starting model, the 2D GL Hamiltonian, takes the 
form 



H = 



d 2 r[e |*(r)| 2 + $|(-iV + 2| e |A)*(r)| 2 (1) 



with the partition function Z = Trexp(— H/T), where 
£o is the coherence length at T — 0, £o = — 1 + T/T c0 , 
\£(r) is the pair-field, s is the film thickness, and the 
Landau gauge A = (0, Hx, 0) will be used hereafter. The 
magnetic screening due to the fluctuation of the gauge 
field A will be neglected based on the familiar reasoning 
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that the effective penetration depth pj A = A 2 /s defined 
in the Meissner state is, in most cases, beyond the system 
size, where A is the London penetration depth. 

Our simulations have been performed by fixing [8j the 
number of vortices N s in the manner commensurate with 
the triangular lattice 



L x L y = 2nr„N, 



H 1 



2 ' 



(2) 



is assumed to be satisfied. Due to this, any gauge- 
invariant quantities such as ^(r)! 2 becomes periodic in 
the perfect vortex lattice. 



The pair-field \& will be expanded in terms of the LL 
eigenfunctions i/Vi.z in the way [B|, Q 



where ru = (2|e|if) 1 / 2 , Lj is the system size in each 

1/2 

direction, and N s is assumed to be an integer. Thus, 
in our simulations, each system size Lj increases with 
decreasing H . Further, the quasi-periodic boundary con- 
dition for the pair- field ^(r) 



f = 



N s -1 



^2 ^2c n ,l1pn,l( T ) 



n=0 I 



(5) 



V(x,y + L y ) 
9(x + L x ,y) 



exp(- 



ir H 2 L x y) 



(3) 
(4) 



Then, eq.© becomes n[bT/s] = H[bT/s]/T, where 



H[bT/s] = - 1 + (2n+ l)M|c„,i| 2 + 



bT 



2sL x Ly 



E 



E E E c nidi c n2,l2 Cn 3,h c n i ,h^li+h+N s (N 1 +N 2 )d2+U 
NiN 2 { ni } {k} 



r H iK x K y r H 



(6) 



Here, fc + = fc^ + ik y , and 
2tt 



2tt 



fcx = -j—m x , k y = —(h-h + N s N-_ 

L X Ly 



(7) 



with integers m x , li, I2, and iVi. Further, the quartic 
term has been rewritten according to the treatment used 
elsewhere and in terms of the expression 



min(p.q) 



n>0 



(p — n)\(q — n)\n\ 



z p - n (-z*) q - n . (8) 



In the regime where the phenomenological GL model is 
applicable, the coefficient b is given, in terms of the GL 
parameter k and the depairing field H C 2(0) at T = 0, by 

As numerical values of material parameters, we use here- 
after those typical of optimally-doped high T c cuprates, 
T c0 = 10 2 [AT], £ = 10[A], and s = 15[A]. Then, 
bT/(L x L y s) is expressed by W(t, h)/(nN s ), where 

W(h,t) = k 2 ■ IQ~ 5 th, (10) 

t = T/T c q, and h = H/H C 2(0). By performing the scale 
transformation 



{Cnj} -> {Cn,l}/y/W(h,t), 



(11) 



eq.© finally becomes 



rlred — 



1 



W(h,t) 



(12) 



with Z = Trexp(— H r ed) which will be used for numerical 
simulations. 

As a thermodynamic evidence of a first order transi- 
tion, we focus on the hysteresis AE of the internal energy, 
which can be defined by [|[ 



AE = 



1 



L X Ly 



(H) dcc - (H)i 



(13) 



where the indices "dec" and "inc" denote cooling and 
warming processes, respectively, and (• • • ) implies the 
thermodynamic average. In the warming process, the 
vortices is assumed to form the triangular lattice at a low 
enough temperature. This condition can be expressed as 

SI 



, »t1/2x 1/2 

Cl M = (-ir (m+1)/2 (^j s lfi s, 



.N 



1/2 (14) 

3 m/2 x ' 



for < m < 2Ns /2 but is zero otherwise, where 0a de- 
notes the Abrikosov factor (= 1.1596) of the perfect tri- 
angular lattice. On the other hand, the initial condition 
c n .i = will be used for the cooling process. 
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III. NUMERICAL RESULTS 

In our simulations, the field dependence has been ex- 
amined at the fixed number of field-induced vortices N a 
and by changing Lj following our previous work in LLL 
Q , because numerical results become much clearer than 
those under fixed Lj and a field- induced change of N s . 
We have performed three simulations with N s = 36 and 
six LLs (0 < n < 5), N s = 64 and the six LLs at only 
h = 0.1 and 0.3, and N s = 36 and seven LLs (0 < n < 6). 
We have not found any ^-dependent essential differ- 
ences in the hysteresis data and snapshots of the vor- 
tex configurations. Hereafter, we will primarily show the 
data resulting from the use of six LLs and N s =36. 

Simulations have been performed by creating Markov 
chains of the coefficients c n ,i according to the Metropolis 
algorithm. MonteCarlo (MC) steps of the range between 
5.0 x 10 5 and 1.0 x 10 6 were used to ensure approach to the 
thermodynamic equilibrium, and, afterward, additional 
5.0 x 10 4 MC steps were used to take a statistical average 
of physical quantities. 

Now, let us explain our numerical results. Figure 1 
expresses the hysteresis, eq . (TIB")) , accompanying the first 
order melting transition at t m {h) = T m (H) /T c q at two 
different magnetic fields, (a) h = 0.3 and (b) h — 0.6, on 
sweeping the temperature in the LLL approximation (up- 
per figures) and in the case with higher LLs with n < 5 
(lower ones). It is found that, compared with the familiar 
LLL results, inclusion of higher LLs depresses t m and re- 
duces the hysteresis around t m . In h = 0.3, these higher 
LL effects are clearly seen, while the LLL approximation 
seems to be valid even quantitatively in h = 0.6. 

Prior to a further discussion on the field dependence of 
the hysteresis, higher LL effects on the h-t phase diagram 
will be explained here. In Fig. 2, the melting transition 




FIG. 1: Numerical t (= T/T c0 ) v.s. AE curves in (a) h = 0.3 
and (b) h = 0.6 taken at the fixed N s = 36. For both, the 
upper figure is the result in the LLL approximation, while the 
lower one is obtained by taking account of the six LLs with 
< n < 5. 
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FIG. 2: Resulting melting transition curve t m {h) (solid curve 
with open circles) from the simulation with six LLs incorpo- 
rated. The dashed curve is the corresponding curve in the 
LLL approximation. The thin dotted curve is the H C 2(T) 
line, while the thick dotted curve is the estimated boundary 
above which the thermally-induced vortex-pairs appear. 



curve (solid curve with open circles) following from the 
present work is shown together with the corresponding 
result (dashed curve) in the LLL approximation. The 
melting temperature at each /i-value has been defined 
as the i-value at which the hysteresis becomes maximal. 
In obtaining these curves, we have used the value k = 
61 for the GL parameter. In LLL approximation, the 
reduced melting temperature t m yields the LLL scaling 
1 — t m — h oc (t m h) 1 / 2 [T3|, and the results in Fig. 2 show 
that 



1 - t r 



-h = c^ 1 



W(t m ,h) 



(15) 



with @ C2 = 0.0989. The r.h.s. of eq. (fl"5| ) is proportional 
to (M) 1 / 2 , and the resulting t-h relation is called the LLL 
scaling [lflj. On the other hand, in low enough fields, the 
relation 



1 - tm - h = 



W(t m ,h) 



(16) 



2nc'h 

is expected to be satisfied [l(| ■ Reflecting the fact that 
the r.h.s. of eg- flB"] ) is /i-independent, the resulting melt- 
ing curve in lower fields is nearly vertical in the t v.s. h 
phase diagram. One can verify that, in Fig. 2, the weak 
field dependence of the solid curve in such low fields is 
a reflection of the field dependence of T c2 (h) rather than 
the neglect of other higher LLs (n > 6). In contrast to the 
LLL approach, however, there is no well-established value 
of c' in eq.flTfi"]). Figure 2 suggests that t m — 0.3 in low 
enough fields, which implies that c' = 0.0025. Namely, 
the relation c' = 0.25c 2 . is approximately satisfied in the 
results shown in Fig. 2. 

Here, we will compare the solid curve in Fig. 2 with the 
corresponding one following from the elastic theory. To 
do this, the results in the elastic theory [l|, [ll| should be 
reviewed. In any approach based on the elastic energy 
of the vortex lattice, the field dependence of the melting 
temperature T m (H) is determined by that of the shear 
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elastic molulus Cqq, which depends on the magnitude of 
the reduced applied field h = H/H C 2(0), where H C 2(T) is 
the depairing field in the orbital limit, even in the low 
field regime where the phase-only model is useful. Then, 
C*66 oc H in lower fields, while ~ {H — H C 2) 2 in higher 
fields [12J. Consequently, T m (H) is iJ-independent in 
low fields, while T m (H) obeys the lowest LL scaling [![. 
Clearly, the low field portion of the solid curve in Fig. 2 is 
comparable with the H- independent T m (H)-cmve in the 
London limit. The above-mentioned < m (i?)-curve in the 
elastic theory can be described by a single expression 



ee r H 
Tco 



(17) 



implying a comparison between the elastic and thermal 
energies by assuming that the mechanism of the vor- 
tex lattice melting is universal and uniquely given irre- 
spective of the magnetic field strength, where the coef- 
ficient a is ^.-independent. The vortex shear modulus 
Cq6 is expressed, in the case of the triangular lattice, 
as H C 2{T) H I (32ir k 2 ) in the low field London regime and 
0.708(# c2 (T)-iJ) 2 /(327rK 2 ) in the high field LLL regime 
[l2|, where the value Pa = 1-1596 for the triangular lat- 
tice was used. Then, we find that the elastic model leads 
to the relation c' = 0.704c 2 , when t m (h) in the above- 



mentioned two regimes is expressed by eqs. ([T5|) and f|X6j> . 
That is, the present t m (H)-iesult shown in Fig. 2 sug- 
gests that the melting temperature in low fields where 
the LLL approximation breaks down is over-estimated in 
the elastic theory. The main origin of this discrepancy 
consists in the presence of higher LL fluctuations, which 
are not included in the elastic theory, in the present nu- 
merical simulations of the GL model. Since the mean 
field solution of the triangular or hexagonal vortex lat- 
tice is described only by the LLs with LL indices of mul- 
tiples of six, reflecting its six-fold orientational symmetry 
[l3j . the shear modulus in this state is also determined 
by those LLs 14|. In other words, the higher LLs with 
1 < n < 5 are not incorporated in the relation (TIT)) . On 
the other hand, the n = 6-LL is not included in our com- 
putation leading to Fig. 2, while the mean-field GL theory 
should reduce to the London theory just by incorporating 
the higher LL modes with LL indices of multiples of six. 
Namely, the nearly vertical t m {h)-\me in lower fields in 
Fig. 2 has been obtained irrespective of reduction to the 
London model due to the lowering of h. In other words, 
the elastic model overlooks crucial fluctuation effects and 
is not sufficient for describing the vortex lattice melting 
as far as choosing the GL model as the starting model is 
valid. 

Returning to Fig.l, one can see that the magnitude 
of the hysteresis in the simulations including higher LLs 
has a different field dependence from that in the LLL ap- 
proximation. The hysteresis curves following only from 
LLL indicate that, as is seen from the upper figures in 
Figs.l and 3, the hysteresis accompanying the transition 



increases with decreasing h. Once the higher LLs are in- 
corporated, however, the hysteresis peak at the transition 
rather decreases with decreasing h, implying a reduction 
of the first order character of the melting transition in 
lower fields. Of course, this field dependence correlates 
with the corresponding /i-dependence of the difference in 
the £ m -value between the LLL curve (dashed curve) and 
the solid one. However, this trend that the first order 
transition is weakened by decreasing the field is inter- 
esting in that it suggests the possibility that, in lower 
fields, the first order melting might be transmuted to the 
two-step continuous melting scenairio To clarify this 
possibility within our GL study, the corresponding hys- 
teresis curves in lower fields are shown in Fig. 3 together 
with the /i-dependence of the peak height of hysteresis at 
t m (Fig. 3(a)). Although the hysteresis reduces with de- 
creasing field even for lower h (= 0.1 and 0.05), the data 
suggest that the reduction of hysteresis saturates with 
vanishing h. This saturation of the peak height does not 
seem to be an artifact due to the limitation to the LLs 
with n < 5: An additional (dashed) hysteresis curve for 
h = 0.1 obtained by including the n = 6 LL further 
is presented in Fig. 3 (a). Although the transition point 
has been slightly shifted to a lower temperasture by the 
n = 6 LL which, as mentioned earlier, affects the elastic 
energy of the mean field vortex lattice, the hysteresis at 
t m is slightly bigger by including the n = 6 LL. There- 
fore, based on the present simulation, we argue that the 
melting transition remains of first order even in low field 
limit, although its first order character is significantly 
weakened with decreasing the field. 

Next, we explain the vortex states at various t- values 
seen in our simulation using the LLs < n < 5. In 
Fig. 4, we have shown snapshots of the distributions of 
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FIG. 3: Numerical t v.s. AE data at (a) h = 0.1 and (b) 
h = 0.05. The upper figure of (a) is the result in the LLL 
approximation, while the lower one of (a) includes the solid 
curve in the case with six LLs (n < 5) and the dashed one in 
the case with seven LLs (n < 6). The figure (c) implies the 
/i-dependence of AE m = AE(t — t m ). Data in (b) and (c) 
result from the use of six LLs with n < 5. 
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t = 0.2[Rp^t = 0.3 



l = 0.2 



3fc S: 



t = 0.6 



FIG. 4: Snapshots of the spatial distribution of |^(r) | (top fig- 
ures) and the gauge-invariant gradient vector, eq. (|18p . (mid- 
dle and bottom ones) at each t indicated in the figures. The 
bottom figure is the zoom of the region specified by a square 
in the middle one at t — 0.6 and shows the presence of one an- 
tivortex (circle) with a clock-wise circular current in contrast 
to others with a counterclock-wise circular current. 



the order parameter amplitude \^>\ (upper figures) and 
the gauge-invariant gradient 



00 



(18) 



of the phase <& of \& (middle and lower ones). Before 
discussing what the figures imply, the relation between 
<F and the vortices in the LLL approximation will be 
reviewed. As is well known, the order parameter within 
the LLL satisfies 151 



Bearing this role of the higher LLs in mind, one finds in 
Fig. 4 where h = 0.1 that, below t = 0.5, the antivor- 
tices do not appear. As is seen in the bottom figure, 
appearance of an antivortex and thus, of one vortex pair 
induced by the thermal fluctuation is verified just above 
t = 0.5. By similarly defining the temperature, at which 
the thermal vortex-pairs begin to appear, at different h- 
values, we obtain the crossover line, shown in Fig. 2 as the 
thick dotted line, which separates the vortex liquid com- 
posed only of the field-induced vortices from the state 
with the thermally-induced antivortices. This result is 
consistent with the picture argued elsewhere [l(| that 
the higher LLs with antivortices included become more 
important at higher temperatures and push the melting 
transition curve down to lower temperatures (see Fig. 2 
(a) in Ref.0). 

Finally, we note that, as is seen at the top of Fig. 4, the 
snapshot taken just at t m (h = 0.1) = 0.3 does not in- 
clude any dislocation. This is strange at least within the 
conventional picture based on the elastic theory that the 
melting transition is driven, more or less, by thermally- 
induced dislocations 4]. From the figure, we feel that 
spatial variations of the amplitude |^| of the pair-field 
assist the nearly harmonic shear elastic modes and that 
their cooperative roles result in the melting of the vortex 
lattice. This view arguing the necessity of the amplitude 
fluctuation does not contradict the conventional wisdom 
that Goldstone modes in an ordered phase cannot be- 
come critical modes for driving a thermal disordering. 
Nevertheless, we have to mention that it is beyond the 
scope of the present work to judge whether such a role 
of the fluctuation of the amplitude |^| is essential or an 
artifact of the use of the GL model which is usually valid 
when |>F I is small. 



IV. SUMMARY AND DISCUSSION 



* L (r) = cxp( -^-g- )f Ns (y 



(19) 



in the present gauge A = Hxy, where the function f m (z) 
is any m-th order polynomial of z. Thus, >Fl has N s 
zero points r„ = (x v ,y v ) and takes the product form 
n„ [y — y v + i(x — x u )] Then, it is staightforward to show 
that the phase $ of "Fl satisfies the topological condition 



N. 



V x V$ = 2tt 2^ S(x - x v )8{y - y v )z, 



(20) 



implying each zero point of is a vortex coordinate. 
Equation (f20| implies that, in all of the N s zero 
points of I^lI corresponds to the field- induced vortices 
and thus that no antivortices can appear. That is, the 
thermally-induced vortex pairs are described not by the 
LLL modes but only by the higher LL modes of <F (loj . 



In the present work, numerical simulations of the GL 
model with not only the order parameter modes in LLL 
but also those in five or six higher LLs have been per- 
formed. It has been found that the nearly vertical melt- 
ing curve, which has been identified so far with the result 
following from the elastic model in the London limit, is 
obtained simply by incorporating the lower five higher 
LLs which do not contribute to the elastic model, and 
that the first order character of the melting transition di- 
minishes with decreasing the field, although the melting 
picture with two continuous transitions [ij is not reached 
even in the low field limit. The weaker first order tran- 
sition in lower fields suggests that the discontinuous na- 
ture of the transition in clean limit is easily lost by a weak 
pinning effect in real systems with quenched disorder and 
thus, explains why the first order melting transition has 
not been reflected, e.g., in transport data in real super- 
conducting films in nonzero field 16 1. 
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In contrast to the present result in the GL model show- 
ing a small but nonzero hysteresis at the melting tran- 
sition in any magnetic field, a simulation work 17j has 
recently been reported indicating a continuous melting 
transition. There, the authors have discretized the re- 
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lation Uj = r 2 H (z x V<5<fr)j, justified [lj, ll8| within the 



linear elastic theory, between the shear displacement u 
and the phase fluctuation <5$ to invoke a starting lattice 
model for their numerical studies. However, their start- 
ing model includes phase slips in the j-direction and the 
direction perpendicular to this on the same footing, sug- 
gesting that thermally-induced vortex pairs coexist with 
the thermally-induced dislocation pairs. That is, the 
model in Ref. [17j is not compatible with the GL model in 
which the vortex pairs never appear close to the melting 
transition (see Fig. 2), and thus, it is not surprising that 
the present results are not consistent with the continuous 
transition [ItJ found in a model which, in our opinion, is 
incompatible with the original GL model. 

On the other hand, we do not definitely conclude at 
this stage the absence of the low field regime in which the 
melting transition is continuous 0] with no thermally- 
induced vortex pairs. More or less, the main drawback 
of the present work is the use of the GL model in ad- 
dressing the low field regime. In such low fields and low 
temperatures, the amplitude |^| is rigid enough to justify 
the phase-only model (with no thermally-induced vortex 
pairs) . On the other hand, the fact that the nearly verti- 
cal melting curve in low fields expected from theoretical 
arguments is obtained simply by taking account of five 
or six higher LLs' modes suggests that the results are 
not significantly affected by the truncation of the num- 
ber of incorporated LLs. Nevertheless, further theoreti- 
cal progress will be necessary to resolve this issue on the 
character of the transition. 



[1] 
[2] 

[3 
[4 

[5. 
[6 

\-\ 

[8 

[9 

[10 
[11 

[12 
[13 
[14 

[15 
[16 



D. S. Fisher, Phys. Rev. B 22, 1190 (1980). 

D. S. Fisher, M. P. A. Fisher, and D. A. Huse, Phys. Rev. 
B 43, 130 (1991). 

T. Nattermann and Scheidl, Adv. Phys. 49, 607 (2000). 
B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 
121 (1978). 

Y. Kato and N. Nagaosa, Phys. Rev. B 48. 7383 (1993). 
J. Hu and A. H. MacDonald, Phys. Rev. B 56, 2788 
(1997). 

P. G. de Gennes, Superconductivity of Metals and Alloys 
(Addison-Wesley, 1989) Ch.3. 

K. Myojin, R. Ikeda, and S. Koikegami, Phys. Rev. B 78, 
014508 (2008). 

N. Hiasa, T. Saiki, and R. Ikeda, Phys. Rev. B 80, 014501 
(2009). 

R. Ikeda, J. Phys. Soc. Jpn. 64, 1683 (1995). 
S. Doniach and B. A. Huberman, Phys. Rev. Lett. 42, 
1169 (1979). 

E. H. Brandt, J. Low Temp. Phys. 26, 735 (1976). 
G. Lasher, Phys. Rev. 140, A523 (1965). 
R. Ikeda, T. Ohmi, and T. Tsuneto, J. Phys. Soc. Jpn. 
59, 1740 (1990). 

Z. Tesanovic, Phys. Rev. B 44, 12635 (1991). 
Alternatively, the absence of a sharp vanishing of the lin- 
ear resistance in 2D superconductors may be attributed 
to the absence of 2D vortex glass phase, as argued else- 
where [R. Ikeda, J. Phys. Soc. Jpn. 65, 3998 (1996)]. 
Effects of the diminishing of the first order nature of 
the melting transition in clean limit would be clarified 
through other probes such as thermodynamic quantities. 
J. Iaconis, R.G. Melko, and A. A. Burkov, Phys. Rev. B 
82, 180504 (2010). 
[18] M. A. Moore, Phys. Rev. B 39, 136 (1989). 



[17] 



